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We discuss a nonlinear model for the relaxation by energy 
redistribution within an isolated, closed system composed of 
non-interacting identical particles with energy levels with 
i — 1, 2, ...,N. The time-dependent occupation probabili- 
ties Pi(t) are assumed to obey the nonlinear rate equations 
r dpi/dt — —pi In pi — a(t)pi — fi(t) api where a(t) and f3(t) 
are functionals of the Pi(t)'s that maintain invariant the mean 
energy E — 53 -=i e, Pi(t) and the normalization condition 

1 = E£U The entr °py s(t) = -k B Y,f =1 Pi(t)^Pi(t) 

is a non-decreasing function of time until the initially nonzero 
occupation probabilities reach a Boltzmann-like canonical dis- 
tribution over the occupied energy eigenstates. Initially zero 
occupation probabilities, instead, remain zero at all times. 
The solutions pi (t) of the rate equations are unique and well- 
defined for arbitrary initial conditions Pi(0) and for all times. 
Existence and uniqueness both forward and backward in time 
allows the reconstruction of the ancestral or primordial low- 
est entropy state. By casting the rate equations not in terms 
of the pi's but of their positive square roots y/pi, they un- 
fold from the assumption that time evolution is at all times 
along the local direction of steepest entropy ascent or, equiv- 
alently, of maximal entropy generation. These rate equations 
have the same mathematical structure and basic features of 
the nonlinear dynamical equation proposed in a series of pa- 
pers ended with G. P. Beretta, Found. Phys. 17, 365 (1987) 
and recently rediscovered in S. Gheorghiu-Svirschevski, Phys. 
Rev. A 63, 022105 and 054102 (2001). Numerical results illus- 
trate the features of the dynamics and the differences with the 
rate equations recently considered for the same problem in M. 
Lemanska and Z. Jaeger, Physica D 170, 72 (2002). We also 
interpret the functionals k B a(t) and k B /3(t) as nonequilibrium 
generalizations of the thermodynamic-equilibrium Massieu 
characteristic function and inverse temperature, respectively. 



I. INTRODUCTION 

Much work has appeared in recent years on the study 
of entropy-generating irreversible nonequilibrium dynam- 
ics. Limited discussions of previous work is found in [1-3] 
and references therein, but no thorough critical review of 
the subject is available, although it would be very helpful 
to provide proper acknowledgement of pioneering work, 
avoid 'rediscoveries' such as in [4] and outline the dif- 
ferent frameworks, motivations, approaches and contro- 
versial aspects. To be sure, recent discussions [2,4-6] on 



possible fundamental tests of standard unitary quantum 
mechanics, related to the existence of 'spontaneous deco- 
herence' at the microscopic level, and on understanding 
and predicting decoherence in important future appli- 
cations [7] involving nanometric devices, fast switching 
times, clock synchronization, superdense coding, quan- 
tum computation, teleportation, quantum cryptography, 
etc. show that the subject of irreversible nonequilibrium 
dynamics is by no means settled. 

It is not the purpose of this paper to attempt such a 
difficult review, nor to address the related fundamental 
issues lurking beneath interpretation (see, e.g., [8-10]). 
Rather we wish to address the model problem recently 
outlined in [1]. 

This model may prove useful to complement vari- 
ous historical and contemporary efforts to extend linear 
Markovian theories of dissipative phenomena and relax- 
ation based on master equations, Lindblad and Langevin 
equations, to the nonlinear and far nonequilibrium do- 
main. For example, spectroscopic studies of the ef- 
fects of vibrational relaxation on line shapes of two- 
level electronic transitions cannot be regularized under 
the Markovian approximations so that various nonlinear 
approaches are being developed and tested [11] in some 
cases at the expense of giving up preservation of (com- 
plete) positivity [12] or hermiticity [13] of the (reduced) 
density operator. 

Again, it is not our purpose here to review the litera- 
ture of these specific potential applications of our model 
dynamics, nor to apply it explicitly to particular exam- 
ples. Rather we wish to focus on illustrating its general 
features (including preservation of positivity and her- 
miticity at all times, even backwards) that make it a 
good candidate (that is, compatible with all reasonable 
requirements imposed by thermodynamic principles [14]) 
of extensions of the traditional linear master equations 
for open system dynamics, capable to include the de- 
scription of nonlinear spontaneous relaxation within the 
system (even if isolated) by energy redistribution between 
the occupied levels. 

We consider an isolated, closed system composed of 
non-interacting identical particles with single-particle en- 
ergy levels ej with i = 1,2, ... ,N where N is assumed 
finite for simplicity and the e^'s are repeated in case of 
degeneracy. We restrict our attention on the class of 
dilute-Boltzmann-gas states in which the particles are 
independently distributed among the N (possibly degen- 
erate) one-particle energy eigenstates. In density oper- 
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ator language, this is tantamount to restricting the at- 
tention on the subset of one-particle density operators 
that are diagonal in the representation which diagonal- 
izes the one-particle Hamiltonian operator. We denote 
by pi the occupation probability of the i-th eigenstate, 
so that the per-particle mean energy, normalization and 
entropy functional are given by the relations 



N 



N 



i=i 



N 



5(p) = -k B ^2p l \np l , 



(1) 



i=l 



where p denotes the vector of p^s, the Boltzmann con- 
stant k B may be used to nondimensionalize S (or we may 
assume for simplicity k B — 1 unit of entropy), and of 
course U(p) = 1 for any normalized distribution p. 

As is well known, for a given value of E, the 
thermodynamic-equilibrium canonical distribution 



pT(e) = 



£Lexp(-/3 se (£) ei ) 

has inverse temperature (3 se (E) = l/k B T(E) and maxi- 
mal entropy S SC {E) = —k B J2iLi PT(E) \npf(E). 

We are interested in studying the dynamics of a 
nonequilibrium distribution obtained, for example, by ex- 
citing some energy eigenstates. As suggested in [1] , a way 
to alter the distribution is to repopulatc (e.g., by selective 
laser heating) or depopulate (in principle, by selective 
cooling or resonance fluorescence) a subset of eigenstates. 
This is described by multiplying each pj by a perturba- 
tion factor fj > (with j = 1, 2, . . . , N) (repopulation 
fj>l, depopulation fj < 1) and then renormalizing, to 
yield the perturbed nonequilibrium distribution 



(2) 



Pj 



fjPf(E) 



ElifiPT(E) 

Of course, in general the perturbed distribution has a dif- 
ferent mean energy E = Eil=i eiPi and different entropy 
S = — k B Ej=i Pi m Pi- However, a proper choice of the 
perturbation factors fi may maintain E = E, in which 
case S < S SC (E) (see Section VII). 

To describe the relaxation towards the new target 
canonical equilibrium distribution p se (E), the dynami- 
cal equation proposed in [1] is, for j = 1, 2, . . . , N, 



(3) 



dt 



= -v [lnpj + a L (p) + b L (p) ej 



(4a) 



where 



«l(p) = 



Mp) 



E, e 4 E j ej In Pj - Ei h Pi E^ef 



^E*e 4 2 -(E^) 

E» ln P» Ej e J"E» e i ln Pi 
^E^?-(e^) 2 



(4c) 



This equation does have the capability of continuously 
rearranging the distribution so that the perturbed distri- 
bution evolves towards the maximal entropy target dis- 
tribution given by Eq. (2) with energy E. However, in the 
far-noncquilibrium region it has the defect to imply the 
unphysical feature that an initially unpopulated eigen- 
state gets populated at an infinite rate. This feature is 
in contrast with a wealth of successful models of physical 
systems in which by limiting our attention to a subset of 
relevant single-particle eigenstates we get good results, 
that are relatively robust with respect to adding to the 
model other less relevant, unpopulated or little populated 
eigenstates. According to Eq. 4, instead, distributions 
where some eigenstates are very little populated would 
survive only for extremely short times. 

The equation of motion that we propose for the time 
evolution of the perturbed distribution is, for j = 
1,2, ... ,N, 



dPj 
dt 

where 



\pj liipj + a(p) P j + /3(p) ejpj] 



(5a) 



Ei dPi Ej ejPj In Pj - E* Pi ^ Pi Ej CjPj 
a(p) = J - ; — J -^— , (5b) 



0(P) = 



E* 4pi - ( Ei e iPi) 

Ej Pi ln Pi Ej e jP] - E» ^Pi ]n P^ 

Y.APi - (jZi^Pi) 



(5c) 



We show in Section II that the apparently slight mod- 
ification with respect to Eq. (4) not only fixes the cited 
defect, while maintaining the relevant overall features of 
conserving energy, normalization, nonncgativity of the 
probabilities and maintaining the entropy generation rate 
nonnegative. It also features existence and uniqueness 
of the solutions of the Cauchy problem for all times, 
— oo < t < +oo, and entails a large class of partially- 
canonical equilibrium distributions that are unstable, as 
well as a single conditionally-stable canonical equilibrium 
distribution for each value of the energy, as required by 
a well-known statement of the second law of thermody- 
namics [15,16]. 

We show in Section II that the structure of Eq. (5) 
is the same as that of the general nonlinear quantum 
equation we discuss in a series of papers written over 
twenty years ago [17-21] in which we develop and propose 
a nonlinear quantum dynamics in an attempt to unite or- 
dinary quantum mechanics and general equilibrium and 
nonequilibrium thermodynamics. As acknowledged also 
in [2,22], the nonlinear quantum dynamical law first pro- 
posed by this author [23] does have very intriguing and 
appealing mathematical features. We must admit how- 
ever that the physical interpretation, motivation and con- 
text of our pioneering scheme is still considered 'adven- 
turous' [24] by most of the physical community, although 
we would prefer to term it 'revolutionary' in the sense 
of Kuhn [25]. For this reason, in this paper we do not 
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pursue such controversial interpretation, but we wish to 
emphasize that — leaving aside its interpretation and fo- 
cusing attention only on its mathematics — our previous 
work represents to our knowledge the first time that the 
steepest-entropy-ascent (or maximal-entropy-generation) 
ansatz has been explicitly formulated and implemented 
in a general dynamical law capable of describing the re- 
laxation of arbitrary nonequilibrium states towards ther- 
modynamic equilibrium. 

In Section III we provide a derivation of Eq. (5) from 
the assumption that the occupation probability distribu- 
tion evolves along the steepest-entropy-ascent trajectory 
in the state space defined in terms not of the pi's but of 
their positive square roots ^/Pi , s. In Section IV we derive 
a fluctuation-dissipation formulation of the equation and 
in Section V a variational formulation. 

In Section VI we discuss a simplest degenerate case 
in which the relaxation equation admits an analytical 
solution, and we compare results with those numerically 
derived from the natural extension of Eq. (4) to such 
case. Finally, in Section VII we show some numerical 
results that illustrate the general features of the proposed 
nonlinear relaxation equation. 



where the non-negativity follows from the well-known 
properties of Gram determinants (see also Section III) . 

Eq. (5) or the equivalent Eq. (6) is well-behaved in the 
sense that the following general features can be readily 
verified (detailed proofs in [4,17]): 

• it conserves the normalization of the distribution 
and the mean energy E along the entire time evo- 
lution; 

• it preserves the non-negativity of each pi ; 

• it maintains the rate of entropy generation non- 
negative at all times; 

• it maintains unoccupied all the initially unoccupied 
eigenstates; in other words, given a distribution p, 
and defining the vector S(p) of Si's such that, for 
each i = 1,2, ... ,N, Si = if pi = or Si = 1 if 
Pi ^ 0, the vector 6 is time invariant; 

• it drives any arbitrary initial distribution p(0) to- 
wards the partially-canonical (or canonical, if Si = 
1 for all pi's) equilibrium distribution, reached as 
t — > oo, 



II. MAIN FEATURES OF THE ASSUMED 
NONLINEAR RELAXATION EQUATION 

By analogy with the dynamical law introduced in 
[17,19,23,26], Eq. (5) may be also written as a ratio of 
determinants in the form 



Pj ln Pj pj e jPj 

Pi ^ pi i Y. e iVi 

J2 e l p l In p l e iPi E e i Pi 



dpj_ 
dt 



1 E e iPi 

J2 e l p l £ ef pi 



(6) 



where |-| = det[-], and r is assumed constant [27] and may 
be used to nondimcnsionalize time (or we may assume 
r = 1 unit of time). 

The resulting rate of entropy generation may be writ- 
ten as a ratio of Gram determinants in the form 



dS_ 
~dt 





Y^Pii^pi) 2 £pilnp 4 £>iPilnpi 




Ep» ln P 


i 1 


J2 e tPi 




J2 eiPi In pi e iPi 


£e?Pi 


T 




l £ e iP* 








E e tP* E e l Pi 





> 



(7) 



Sj exp(- 



Eti^^M-P pe (E,S)ei) 



(8) 



where, of course, S = S(p(0)) and the value of (3 pe is 
determined by the initial state through the relation 
E^Iie, pT( E ,8) = E = E(p(0)). Distributions 
(8) are those for which dp/dt = 0, i.e., that satisfy 
the equilibrium condition Pi In pi = —api — fieiPi 
for all'i's and some scalars a and (3. 

Moreover, Eq. (6) is well-behaved not only in forward 
time but also in backward time, consistently with the 
strongest form of the principle of causality, by which fu- 
ture states of a strictly isolated system should unfold dc- 
tcrministically from initial states along smooth unique 
trajectories in state domain defined for all times (fu- 
ture as well as past). Indeed, for any given arbitrary 
'initial' distribution p(0) we can follow the unique tra- 
jectory p(t) for — oo < t < +oo. In forward time the 
target distributions of all trajectories are given by Eq. 
(8), p(+oo) = pP°(£(p(0)),<5(p(0))). The backward- 
time earliest (or 'primordial') lowest-entropy distribution 
p(— oo) is also uniquely identified by the given initial dis- 
tribution p(0) through Eq. (6), but it is harder to char- 
acterize analytically in general. Depending on the given 
p(0), the redistribution among energy eigenstates may af- 
fect some of the occupation numbers in a non-monotonic 
way. In the limit as t — > — oo, however, all dpj/dfs [the 
rhs of each of Eqs. (6)] become sign-definite; for example, 
they may become all positive except for a particular one 
which tends to p^-(— oo) = 1, so that all others tend to 
zero, p. oo) = [this can happen only if the mean 

energy E(p(0)) is exactly equal to the fc-th energy level, 
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i.e., only if E(p(0)) = e^], or they may all tend to zero 
except for two particular Pj's, say pj and pj; which tend 
to finite values, clearly with pj + p-^ = 1 (examples in 
Section VII). 

Because the model equation maintains the rate of en- 
tropy generation non-negative, the entropy functional 
S [Eq. (1)] is an 5-function [16] and, therefore, every 
thermal-like canonical equilibrium distribution (8) is <5- 
-E-conditionaHy stable, that is, stable with respect to per- 
turbations that do not alter the mean value E of the 
energy and the set of unoccupied energy eigenstates (de- 
scribed by the zeroes in vector 8). These distributions 
constitute the 'target' highest-entropy states compatible 
with the mean value of the energy and the invariant sub- 
set of unoccupied eigenstates. These distributions, how- 
ever, are not ^-conditionally stable, that is, stable with 
respect to all perturbations that do not alter the mean 
value E. Indeed, starting from a distribution (8), a per- 
turbation that changes a zero probability to an infinites- 
imal value, makes the perturbed distribution proceed in 
time by amplifying that probability until a new, differ- 
ent and higher-entropy, target canonical distribution is 
reached. For a given mean energy E, the only canoni- 
cal distribution that is ^-conditionally stable is the one 
for which all energy eigenstates are occupied, i.e., the 
maximal-entropy canonical distribution (2). 

By interpreting the entropy S as a measure of how 
'well' the energy is distributed among the available en- 
ergy eigenstates, the proposed nonlinear dynamics de- 
scribes a spontaneous internal redistribution of the en- 
ergy along the path of maximal entropy increase lead- 
ing towards an 'optimally' distributed (highest-entropy) 
state compatible with the condition of maintaining un- 
occupied the initially unoccupied energy eigenstates. 



III. CONSTRUCTION OF THE EQUATION 
FROM THE STEEPEST-ENTROPY-ASCENT 
ANSATZ 

In this section, we provide a brief derivation of Eq. 
(6) from the assumption that the occupation probabil- 
ity distribution evolves along the steepest-entropy-ascent 
trajectory in the proper state space. We also discuss an 
important degenerate case. 

For the purpose of this derivation, instead of working 
with the vector p of the occupation probabilities, we work 
in terms of their positive square roots [28] , yi — ^Jpl and 
the corresponding vector y. We rewrite the mean energy, 
normalization and entropy functionals as 



N 



N 



4, 



i=i i=i 

N 

s(y) = -k B J2y?^vt (9) 



where, of course, U(y) = 1 for any y. The i-th com- 
ponent of the gradients of these functionals are, respec- 
tively, 

e'i = 2 , u' i = 2y i , s • = -2fc B (^ In y} + y { ) (10) 

and, therefore, the time-rate-of-change functionals may 
be written as 



E = (y,e') , 
S = (y,s') , 



£/ = 2(y,y) , 



(11) 



where (•, •) denotes the scalar product of two vectors [e.g., 
the normalization condition U(y) = 1 may be rewritten 
as (y>y) = 1]> an d the energy and entropy gradient vec- 
tors e' and s' are defined by the components in (10), 
while we choose to substitute immediately the obvious 
relation u' = 2y. 

In order to maintain (y, y) = and (y, e') = the vec- 
tor y must be orthogonal to the linear manifold spanned 
by y and e'. 

For unconstrained maximal entropy generation, y 
would be in the direction of the gradient s' of the en- 
tropy functional S(y); in such case, however, because s' 
is almost never orthogonal to the ye' manifold, in gen- 
eral U(y) and E(y) would not remain time invariant. 
Instead, we assume constrained — constant E(y) and 
U(y) — maximal entropy generation. We obtain it by 
taking y in the direction of the component of s' orthog- 
onal to the ye' manifold. Denoting such component by 
s '±ye' we therefore assume 



4/c B r(y) 



±ye' j 



(12) 



where r(y) may be any positive definite functional of y 
with dimensions of time, that determines the time rate 
at which y evolves along the path of constrained steepest 
entropy ascent. For simplicity, and for the purpose of 
comparison with [1], we assume r a positive constant as 
done in our first proposal of this equation of motion in 
[17,20,21,23]. 

Using the well-known theory of Gram determinants, 
we can write an explicit expression for s'^ ye <. If y and 
e' are linearly independent, we have 



-Lye' — 



(s',y) (y,y) (e'.y) 
(s',e') (y,e') (e',e') 



(y,y) (e',y) 

(y,e') (e'.e') 



(13a) 



If instead y and e' are linearly dependent, i.e., if e' — 2e y 
for some scalar e, the expression is 
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-Lye' 



(s',y) (y,y) 



/(y,y) = s'-(s',y)y , 



where we use (y, y) 
that 



E 
S 



(y,e') 
(y,s') 



(13b) 

= 1. In either case, we readily verify 

= 0, f/ = 2(y,y)=0, 
4rfc B (y,y), (14) 



from which we see that the rate of entropy generation is 
related to the norm of y and is positive definite. 
Combining Eqs. (12) and (13a) we find 



4fc B r y = s' - o(y) y - b(y) e' 



where 



a(y) = 
b(y) = 



(s',y)(e',e')-(s',e')(e',y) 
(y,y)(e',e')-(y,e')(e',y) 
(s',e')(y,y)-(s',y)(y,e') 
(y,y)(e',e')-(y,e')(e',y) 



(15a) 

(15b) 
(15c) 



and, setting back pi — yf and pi = lyijji we readily 
obtain Eq. (5) and the identities a(p) = 1 + a(y)/2fc B 
and f3(p) = b(y)/k B . 

Similarly, combining Eqs. (12) and (13b) we find 



4/c B ry = s' - (s',y)y 



(16a) 



and, therefore, in the degenerate case of y and e' lin- 
early dependent, the relaxation equations are, for j = 
1,2,..., TV, 



or, equivalently, 

d Pj = 
dt 



Pj ln Pj 



Pj 
1 



(16b) 



(16c) 



In this case, the rate of entropy generation may be writ- 
ten as 



dS_ 
~dt 



Y.Pii^Pi) 2 x>;inp; 

J2 Pi In Pi 1 



>o, 



(17) 



where the nonegativity follows from the well-known prop- 
erties of Gram determinants. 

Equation (16c) substitutes (6) when e' = 2ey for some 
scalar e or, equivalently, when e-ip-i — ep t for every i, 
that is, when the populated eigenstates all correspond 
to the same energy level. If satisfied at one instant in 
time this condition is satisfied at all times, both forward 
and backward in time. It follows that in such degenerate 
cases the entire time evolution is governed by Eq. (16b). 



In the general nondegenerate cases, i.e., when at one time 
(and, hence, at all times) eiPi ^ Epi for two or more i's, 
where E is the mean energy, the time evolution is entirely 
governed by Eq. (5) [or, equivalently, (6) or (15)]. This 
also implies that the denominators of a(p) and /3(p) [or, 
equivalently, of a(y) and b(y)] remain positive definite 
at all times and, hence, the entire time evolution is well- 
defined. 

As stated above, the feature that unpopulated eigen- 
states remain unpopulated is extremely important as it 
is compatible, for example, with the widely accepted and 
successful possibility to describe real systems by means 
of simplified models with a limited number of relevant 
energy eigenstates. 

The same feature does not hold for equation (4), be- 
cause it implies that an initially unpopulated eigenstate 
gets populated at an infinite rate. For the same reason, 
Eq. (4) docs not allow tracing the time evolution back- 
ward in time beyond the instant when the first eigenstate 
becomes unpopulated, for at earlier times the condition 
Pi > is not satisfied. 



IV. FLUCTUATION-DISSIPATION 
FORMULATION 

It is noteworthy that Eq. (6) admits a general 
fluctuation-dissipation formulation and interpretation. 
To see this, we introduce the energy and entropy fluc- 
tuation functionals as follows 



(A£A£)(p)=£. Vl [e t -E(p)} 2 

= Pi e 1 - (Ei Pi e * 

(y,y) (e',y) 

(y,e') (e',e') 
(ASAS)(p) = Y,i P* Mb Inpi - S(p)} 2 

= fc B E 2 Pi ( ln ^) 2 - & b(e, Pi ln ^) 2 

(AEAS) (p) = ^ pi [<* - E(p)} [-k B \n Pl - S(p)} 
and rewrite functionals /3(p) and a(p) as 
1 (AEAS)(p) 



(18) 



(19) 
(20) 



0(P) 



k B (AEAE)(p) 



a(p) = ^P--P(p)E(p) . 



(21) 
(22) 



At the thermodynamic equilibrium distribution with 
energy E, p sc (E), we have 

k B p(p sc (E)) = k B (3 sc (E) = 1/T(E) , (23) 
k B a(p sc (E)) = k B a sc (E) = S SC (E) - E/T(E) , (24) 

and in Eq. (24) we recognize the thermodynamic- 
equilibrium Massieu characteristic function [15] 
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M sc = S - E/T 



(25) 



It is therefore natural, in this framework, to adopt the 
following generalization of the Massieu function to arbi- 
trary nonequilibrium distributions 



M(p) - k B o(p) = 5(p) - k B f3(p)E(p) , 



(26) 



with /3(p) given by Eq. (21). The corresponding fluctua- 
tions functional is 

(AM AM)(p) = ^ [-fe lnp 4 - fc B /3(p)e 4 - M(p)] 2 . 

(27) 

We can readily verify that [Eq. (7)] may be rewritten as 
dS 1 



dt k B r 



(AM AM) 



(28) 



and, therefore, the rate of entropy generation is di- 
rectly proportional to the fluctuations of our general- 
ized Massieu function [29] . Such fluctuations are related 
to entropy and energy fluctuations through functional 
&b/?(p) as follows 

(AMAM)(p) - (ASAS)(p) - klp(vf(AEAE)(p) , 

(29) 

and become zero at every canonical thermodynamic- 
equilibrium distribution p se (E), Eq. (2), and at every 
partially-canonical equilibrium distribution p pc (E 7 S), 
Eq. (8), as well. 

It is noteworthy that the functional fc B /3(p), which 
is well-defined by Eq. (21) only for distributions with 
(AEAE) 0, may be interpreted in this framework as 
a natural generalization to nonequilibrium of the inverse 
temperature, at least inasfar as for t — > +oo it tends 
to the thermodynamic-equilibrium inverse temperature 
k B (3 sc of distribution (2) or the partial equilibrium in- 
verse temperature k B (3 pc of distribution (8). 

The special case of distributions with (AEAE) = 
happens if and only if e' = 2ey for some scalar e (see Sec- 
tion III). In such special degenerate case, (AEAE) re- 
mains zero along the entire time evolution, which is given 
by Eq. (16b), and the role of the Massieu function is taken 
up by the entropy 5*, for both equilibrium and nonequi- 
librium distributions. The fluctuation-dissipation form 
of the rate of entropy generation [Eq. (17)] becomes the 
following 



dS_ 
~dt 



(ASAS) , 



(30) 



and in this special degenerate case the canonical and 
partially-canonical equilibrium distributions all have 
(ASAS) = 0, for they consist of N s = (S, 6) proba- 
bilities pi all equal to 1/N& and of N — Ns all equal to 
zero. 



As regards the fluctuation-dissipation relations, the 
various well formulated arguments, derivations and in- 
terpretations discussed for Eq. (4) by Englman in the 
Appendix of Ref. [1] and based on the steepest-entropy- 
ascent ansatz - first introduced in quantum thermo- 
dynamics by the present author [19] - apply with mi- 
nor modifications also for our better-behaved dynamical 
equation Eq. (5). In addition, we prove in [26] that Eq. 
(5) implies a generalized Onsager reciprocity theorem. 



V. VARIATIONAL FORMULATION 

In terms of the yi — ^/pl notation, we can derive our 
equation of motion also as a result of the following equiva- 
lent variational formulation (along the lines recently pro- 
posed in [4]) 

max S = (y, s') subject to E = (y, e') = 0, 
y 

£7=(y,y)=0, and (y, y) = £(y) , (31) 

where the last constraint implies that we maximize the 
entropy generation rate only with respect to the 'direc- 
tion' of y, i.e., at every given y we select the maximizing 
y among a subset of vectors that share the same (but 
otherwise arbitrary) norm £ (y) . For y and e' linearly in- 
dependent, using the standard method, we associate the 
Lagrange multipliers a, b and 4fc B r with the constraints, 
and from Eq. (14) and the necessary Euler-Lagrange con- 
ditions 

A [(y, s') - a (y, y) - b (y, e') - 4fc B r (y, y)] = , (32) 

we readily obtain Eq. (15) as well as, upon substitution 
into the constraints, the multipliers given by Eqs. (15b) 
and (15c), and the square norm of y, 



ay) 



s 



1 r(s',y,e') 



4fc B r 


16fc2r 2 r(y,e') 






(s',s') (y,s') (e', 






(s',y) (y,y) (e', 


y) 


1 


(s',e') (y,e') (e', 


e') 


16fcir 2 


(y,y) (e',y) 






(y,e') (e',e') 





(33) 



where T denotes the Gram determinant of the argument 
vectors. 

Similarly, for the degenerate cases with y and e' lin- 
early dependent, the normalization and constant energy 
conditions collapse into a unique constraint with which 
we associate the Lagrange multiplier c, and by the same 
standard procedure we obtain Eq. (16a) and, upon sub- 
stitution into the constraints, the multiplier c = (s',y), 
and the square norm of y, 
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£(y) = 



s 



1 



4fc B T 16fc2T 2 " r ^' y ^ 

(s',s') (y,s') 



1 



16fc2 T 2 



(s',y) (y,y) 



(34) 



VI. SIMPLEST CASE: TWO-LEVEL PARTICLES 
WITH DEGENERATE EIGENSTATES 

The simplest mathematical form of the model equation 
that derives from the equation of motion proposed in the 
previous sections, is obtained when we have an isolated, 
closed gas composed of non- interacting identical two-level 
particles with degenerate energy levels, such as electonic 
spins in the absence of an applied magnetic field. Then 
N = 2, both levels have energy e\ = e2 — e, the occu- 
pation probabilities of the two corresponding eigenstates 
are p\ = 1 — p and p2 = p, respectively, and the model 
equation (16c) for redistribution among the two eigen- 
states becomes 



dp f 1-p 

— = p 1 - p In 

at p 



(35) 



where we set r = 1. The rate of entropy generation 
(k B = 1) is 



By analogy, the extension to this degenerate case of the 
model equation proposed in [1] is dpj/dt = v(— hxpj+ai) 

with 2cll = St=i m P* = m P + m (l ~ P) tnat i s (setting 
v = 1/2 and adding, for clarity, the subscript L) 



dpL 
dt 



4 PL 



PL 



which yields the entropy generation rate 

2 



S, = i I In 



I -PL 
PL 



(40) 



(41) 



Figure 1 shows a comparison between the time depen- 
dence (38) implied by our rate equation (35) and that ob- 
tained by numerical solution (by a standard Runge-Kutta 
method) of the rate equation (40), for bothp(— oo) = = 
Pl(0) and p(— oo) = 1 = pz,(0). For the purpose of com- 
parison, time t = is selected where pl(0) = or 1 and 
the initial state p(0) is selected so as to emphasize that 
in the limit as t — » +oo we have p(t) sa ph{t)- In fact, we 
readily verify from both Eqs. (35) and (40) that the two 
time dependences have the same asymptotic behavior as 
they approach the equilibrium distribution, that is, 



1 



dp 

dt~2~ P 



and 



(42) 



S=p(l-p) In 



1-p 



(36) 



Not only equation (35) is well-behaved at all times (ex- 
istence and uniqueness of the solution for any initial p(0) 
with < p(Q) < 1), but it can also be integrated to yield 



t = 



p(t) 



dp 



In 



p(o) p (1-p) In 



1-p 



= In 



In 



1-P(0) 

Kg) 

p(t) 



(37) 



or, equivalcntly. 
P(t) - 



i-p(Q) \ expi ~ t/T) 

P(0) J 



1 



1 1 



■ tanh 



-exp(-t/r)hi 

2 PV 7 ' p(Q) 



2 2 

from which we readily hnd p(oo) = 1/2 and 



(38a) 



(38b) 



p(-oo) = 



1 



1-P(0) 
p(0) 



30 



forp(0) < 1/2 

1 forp(0) > 1/2 



(39) 




FIG. 1. Comparison of the time dependences p(t) and pl (i) 
respectively implied by the rate equation (35) that we propose 
and the rate equation (40) discussed in [1]. 



This simplest case brings out the evident different be- 
havior at early times and the unphysical feature of the 
solution of Eq. (40) at pl = where the repopulation 
rate is infinite, implying that no unpopulated eigenstate 
can survive unpopulated. Instead, our Eq. (35) main- 
tains unpopulated any initially unpopulated eigenstate, 
and it also maintains relatively little populated an ini- 
tially little populated eigenstate for a lapse of time that 
is quantified by Eq. (37) and depends on how close the 
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initial value p(0) is to zero. For example, the time re- 
quired to take p(0) = 1CT 2 " to p(t) = 1CT 2 is t « Inn 
(that is, t « r Inn). 




0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 



occupation probability, p, Pl 

FIG. 2. Comparison between the entropy generation rate 
S versus p as given by Eq. (36) for fc B = 1 and r = 1 and Sl 
versus pl as given by Eq. (41) for k B = 1 and w = 1/2. 

Figure 2 shows a plot of the entropy generation rate <S* 
versus p obtained from Eq. (36) compared with Sl versus 
Pl as obtained from Eq. (41), where again the essential 
differences for small values of p and 1—p are singled out. 



For a four-level nondegenerate system, Figure 3 rep- 
resents on the diagram the families of possible canon- 
ical (2) and partially-canonical (8) equilibrium distri- 
butions which in our dynamics are the only ones with 
zero entropy generation rate. We recall that the slope of 
these curves is related to the parameter [3 pe (E, 8) because 
dS pe (E,S)/dE\ s = k B pv c (E,8), which for the canoni- 
cal distribution (all <5,'s equal to unity) is dS(E)/dE = 
k B P(E) = 1/T(E). 




a 0.2 0.4 0.6 0.8 1 1.2 1.4 

Entropy, S/k B 



FIG. 3. Representation on an energy versus entropy dia- 
gram (for N = 4 and nondegenerate eigenstates with energies 
e = [0, 1/3, 2/3, 1]) of the families of possible canonical and 
partially-canonical equilibrium distributions which in our dy- 
namics are the only ones with zero entropy generation rate. 
For example, a horizontal line at E — 0.4 intersects seven 
different families of partially canonical states. 



VII. NUMERICAL RESULTS 

The energy versus entropy diagram introduced by 
Gibbs represents the intersection with the E-S plane 
of the E-S-V-n surface representing the stable ther- 
modynamic equilibrium states of a system, assuming 
that the energy eigenvalues depend on the volume V 
and the amounts of constituents n, so that the sur- 
face is represented by the so-called fundamental relation 
S = S(E, {ej(V, n)}). In [15] the use of such diagram 
has been extended to include the projection onto the 
E-S plane of all other states, i.e., not only the stable 
equilibrium states but also the non-equilibrium and the 
non-stable equilibrium states, with given fixed values of 
V and n and, therefore, a given fixed set of energy eigen- 
values. On such diagram, therefore, one point represents 
in general a multitude of distributions, except at every 
point of maximal entropy for each given value of E (V 
and n are fixed) which corresponds to a unique canoni- 
cal distribution (2), i.e., a unique stable thermodynamic 
equilibrium state. 



The number of possible distributions that share a given 
pair of values of E and S is in general an (N — 3)-fold 
infinity except at maximal entropy for each value of E, 
where the distribution is unique, and at few other notable 
exceptions such as at minimal entropy for each given 
E where the distribution may be unique or sometimes 
many-fold. For all possible distributions represented by 
a given point on the E-S diagram, we may evaluate the 
rate of entropy generation dS/dt according to Eq. (6) and 
select the highest value, that we denote by S max (E,S). 
The result of this numerical computation is sketched in 
Figure 4 where the iso-^max contour curves are plotted 
on the entire allowed domain on the energy versus en- 
tropy diagram (of course, under the restriction to the 
subset of states specified in the Introduction). 
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0.2 0.4 0.6 0.8 1 1.2 1.4 

Entropy, S/ks 



FIG. 4. Representation on an energy versus entropy dia- 
gram (for N = 4 and nondegenerate eigenstates with ener- 
gies e = [0, 1/3, 2/3, 1]) of the iso-Smax contour curves where 
5max represents at each point in the diagram the highest value 
of the rate of entropy generation dS/dt according to Eq. (6) 
among all the possible distributions represented by that point. 

The next Figures show typical time dependences of 
the occupation probabilities that result from the numer- 
ical integration (by means of a standard Runge-Kutta 
algorithm) of Eq. (6) both in forward and backward 
time. All trajectories in these Figures refer to a sys- 
tem with N = 4 and nondegenerate eigenstates with 
e = [0,1/3,2/3,1], and all have the same mean energy 
E = 2/5; they all tend, of course, to the canonical dis- 
tribution p so (2/5) = [0.3474,0.2722,0.2133,0.1671] that 
has inverse temperature /3 se (2/5) = 0.7321. They are 
obtained by assuming for all cases an initial distribu- 
tion p(0) obtained by perturbing the canonical distribu- 
tion p se (E) [Eq. (2)] according to Eq. (3) with the en- 
ergy preserving perturbing factors defined as follows, for 
j = l,2,. ..,7V, 

ft = 1 - A + A with < A < 1 , (43) 

where A is otherwise arbitrary, and also 8 is arbitrarily 
chosen among the possible vectors of 0's and l's com- 
patible with the given value of E and form (8) of the 
distribution p pc (E,8) (see Figure 3), where (3 pe (E,8) is 
computed by solving the relation JT Pi^i^, 8) — E. For 
all subsequent Figures we use A = 0.9. 

Figure 5 shows the time dependence of the occupa- 
tion probabilities that results under the assumptions just 
cited using E = 2/5, A = 0.9 and 8 = [1, 1,0, 1] in Eq. 
(43) and subsequently substituting in Eqs. (3), that is, 

p(0) = ApP°(£,<*) + (l-A)p sc (£) . (44) 




dimensionless time, tjr 



FIG. 5. Top: typical time dependences of the occupation 
probabilities that result from the numerical integration of 
Eq. (6) both forward and backward in time, for N = 4, 
e = [0,1/3,2/3,1], energy E = 2/5, initial state at t = 
from Eq. (44) with A = 0.9 and S = [1,1,0,1]. The dots 
on the right represent the maximal entropy distribution; the 
dots at the left represent the lowest-entropy or 'primordial' 
distribution; the dots in the middle represent the p pe (E,8) 
distribution used in Eq. (44) to select the t = state, plotted 
at the instant in time when the entropy of the time-varying 
trajectory is equal to the entropy of the p pc (i5, S) distribution. 
Bottom: the corresponding time dependence of the entropy 
(left axis) and the entropy generation rate (right axis) . 

It is noteworthy that when the trajectory gets very 
close to the partially-canonical unstable-equilibrium dis- 
tribution pP°(£; = 2/5,8 = [1,1,0,1]) the entropy sur- 
face presents a local 'plateaux' and the entropy gener- 
ation rate drops almost to zero, but shortly after the 
trajectory bends in a direction of steeper slope that 
drives the generation up again until the canonical distri- 
bution p sc {E) = [0.3474,0.2722,0.2133,0.1671] is finally 
approached, with inverse temperature /3 so (2/5) = 0.7321. 
Of course, the entropy is a monotonically increasing func- 
tion of time along the entire trajectory. 

Figure 6 shows the same trajectory as well as six 
other trajectories, but instead of plotting the time de- 
pendence of the occupation probabilities we plot them 
against entropy. The initial (time t = 0) distribution 
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used to obtain these seven sample trajectories are ob- 
tained from Eq. (44) with E = 2/5, A = 0.9 and each 
of the seven partially canonical states corresponding to 
the given value of the energy. These seven states are eas- 
ily identified on the E-S diagram in Figure 3 by draw- 
ing a horizontal line at E = 0.4. For the first, third, 
and sixth trajectories we use the p pe (E, S) states with 
S = [1,0,1,0], S = [1,0,0,1] and 6 = [0,1,0,1], respec- 
tively, which [as apparent from the subsequent Figure 7] 
are lowest-entropy boundary points of the entropy sur- 
face for the given energy, and turn out to be also the 
'primordial' states of the corresponding trajectories. For 
the remaining trajectories we use the p pe (E, S) states 
with S = [1,1,1,0], S = [1,1,0,1], S = [1,0,1,1], and 
S = [0,1,1,1], respectively. These too are boundary 
points of the entropy surface, but they correspond to 
partial maxima (over the subset of distributions with one 
unoccupied eigenstate as specified by the corresponding 
zero element of 6). It is seen that these partial maxima 
affect the trajectories passing nearby by acting as par- 
tial attractors especially in the initial phase of the time 
evolution. 

Figure 7 is a more elaborate representation of the same 
seven trajectories. They are shown four times from dif- 
ferent perspectives on the backgorund of contour plots 
of the entropy surface, for four pairs of occupation prob- 
abilities. Indeed, for N = 4 and fixed energy E, the 
number of independent occupation probabilities is two. 
Thus for four pairs of probabilities (pi~P2, P2~P3,P3~P4, 
P4-P1), we draw the contour plot of the entropy sur- 
face over the entire domain of allowed values (which of 
course are contained in a triangular region of the first 
quadrant), and over this plot we draw the seven trajec- 
tories (and the seven partially canonical states used to 
choose them). To save space, we then rotate each of 
the four graphs (respectively by 45, 135, 225, 315 de- 
grees) and combine them on the same graph in Figure 
7. The figure visualizes clearly that the trajectories in- 
deed follow paths of locally-steepest-entropy-ascent and 
unfold smoothly also backward in time to the 'primor- 
dial' states. We also note that these lowest-entropy states 
exhibit a singular behavior in that, for example, state 
[2/5,0,3/5,0] is the primordial state for two entirely dif- 
ferent trajectories, state [3/5,0,0,2/5] for three other, and 
state [0,9/10,0,1/10] for the remaining two. Moreover, 
the partially canonical states appear as partial attractors 
of trajectories passing nearby, as seen quite clearly for the 
second, fourth and fifth trajectory of Figure 6, which are 
partially attracted by the partially canonical states with 
S = [1,1,1,0], S = [1,1,0,1] and S = [1,0,1,1], respec- 
tively. 




0.6 0.7 0.8 0.9 1 1.1 1.2 1.3 1.4 
entropy, S(t)/k B 




0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 1.1 1.2 1.3 1.4 
entropy, S(t)/k B 



FIG. 6. Plots of Pi(t) versus S(t) for seven sample time de- 
pendences of the occupation probabilities that result from the 
numerical integration of Eq. (6) both forward and backward 
in time, for different initial distributions. 
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FIG. 7. Each rotated quadrant of the graph represents, for 
the corresponding pair of occupation probabilities, a plot of 
the seven trajectories shown in Figure 6 drawn over contour 
plots of the entropy surface. 



VIII. CONCLUSIONS 

The model we propose for the description of the time 
evolution of the occupation probabilities of a perturbed, 
isolated, physical system with single-particle eigenstates 
with energies ej for i = 1, 2, . . . , N, is in good agree- 
ment with general thermodynamic requirements such 
as energy conservation, conservation of normalization 
and non-negativity of the probabilities, entropy nondc- 
crease, ^-conditional stability of the maximal-entropy 
canonical equilibrium states, ^-conditional non-stability 
of each non-maximal-entropy partially-canonical equilib- 
rium states, and existence and uniqueness of solutions for 
all initial perturbed distributions, both in forward and 
backward time. As in our previous work [17,19,20,23,26], 
the proposed rate equations implement the fundamental 
ansatz that nonequilibrium time dependence follows the 
path of steepest-entropy-ascent (or, using the terminol- 
ogy adopted in [2,4], maximal entropy generation). 

The model can be readily generalized to include addi- 
tional constraints and therefore adapted to other physical 
and nonphysical (e.g., information theoretical, biological) 
problems that obey the same maximal entropy formalism 
and the maximal entropy generation rate ansatz. Us- 
ing the formalism developed in Section III it can even 
be readily generalized to different entropy functional or 
nonlinear objective functionals that may be relevant in 
many other contexts that share with the present the basic 
mathematical framework. 
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